# include "scene.h"
# include <stdio.h>
# include <stdlib.h>
# include <errno.h>
# include <string.h>
// # include <math.h>

extern char **parse_str(const char *str, char c, unsigned short start, unsigned short end, unsigned short *nsubs);
extern double str2num(const char *str, short nsubs);

extern void init_cond(struct grid *g)
  {
   unsigned short i, j, k, isec, tn, nr;
#  ifdef READ_RECOV_FILE
   FILE *rptr;
   char *str, **s, recovf[MAX_FNAME_LEN];
   unsigned short n, m;
   int _errno;
#  else
   double ***u, ***v, ***w, ***T;
#  endif

#ifndef READ_RECOV_FILE
   g->t = t_INIT;
   g->dt = DELTAt;

/* Initializing the SECTIONs 0 ~ 2. */
   for(isec = 0; isec < Nsec; isec++)
   for(i = 0; i < g->nr[isec]+2*NrG; i++) for(j = 0; j < g->nth+2*NthG; j++) for(k = 0; k < g->nz[isec]+2*NzG; k++)
        {
         g->u[isec][i][j][k] = /*sqrt(i+j*k);*/(double)0;
         g->v[isec][i][j][k] = /*sqrt(sqrt(i*j+k));*/(double)0;
         g->w[isec][i][j][k] = /*sqrt(sqrt(sqrt(i*k+j)));*/(double)0;
         g->p[isec][i][j][k] = P_REF - g->density0*GRAVITY_DOWNWARD*0.5*(g->z[isec][k+1]+g->z[isec][k]);
#        ifdef SOLVE_ENERGY
         g->T[isec][i][j][k] = T_INIT;//+sqrt(i+j+k);
#        endif
        }

   u = g->u[2]; v = g->v[2]; w = g->w[2]; nr = g->nr[2];
#  ifdef SOLVE_ENERGY
   T = g->T[2];
#  endif
/* Data at outlet-boudary of the SECTION-2 for PREVIOUS time-steps. These are used
   for the applciation of Orlansky BC to determine phase velocities at outlet. */
   for(i = 0; i < nr; i++) for(j = 0; j < g->nth; j++) for(k = 0; k < zCells_BACKUP; k++) for(tn = 1; tn < tn_BACKUP; tn++)
       {
     /* tn = 0: t = t; It is filled later in "orlansky_outlet()" using current time-step data.
        tn = 1: t = t - dt;
        tn = 2: t = t - 2*dt; etc.

        Initially, quantities at both the previous time-steps are kept same as at the current time-step. */

        g->outBC.u[tn][i][j][k] = u[NrG+i][NthG+j][NzG+k-1];
        g->outBC.v[tn][i][j][k] = v[NrG+i][NthG+j][NzG+k-1];
        g->outBC.w[tn][i][j][k] = w[NrG+i][NthG+j][NzG+k-1];
#       ifdef SOLVE_ENERGY
        g->outBC.T[tn][i][j][k] = T[NrG+i][NthG+j][NzG+k-1];
#       endif
       }

#else

   sprintf(recovf, "%02d.%s.dat", RUN_CNT-1, RECOV_FILE);

   if((rptr = (FILE *)fopen(recovf, "r")) == (FILE *)0)
       {
        _errno = errno;
        fprintf(g->lptr, "ERROR::  In 'init_cond()': Unable to open the recovery-file '%s' to read.\n\tSystem errno = %d.\n", recovf, _errno);
        exit(-1);
       }

   str = (char *)malloc((MAX_RECOVF_LINELEN+1)*sizeof(char));

/* Skip first two line what are comments. */
   for(m = 0; m < 2; m++) while(fgetc(rptr) != '\n');

/* Reading "t" and "dt". */
   n = 0; while((str[n] = fgetc(rptr)) != '\n') n++; str[n] = '\0';
   s = parse_str(str, '\t', 1, 2, &n);
   if(n != 2){
              fprintf(g->lptr, "ERROR: In 'init_cond()': Error while reading time and time-step from the recovery file '%s'\n", recovf);
              exit(-1);
             }
   g->t = str2num(s[0], 3); free(s[0]);
   g->dt = str2num(s[1], 3); free(s[1]);
   free(s);

/* Load initial condition data for internal parts and ghost-cells of the SECTION-0~2. */
   for(isec = 0; isec < Nsec; isec++)
     {
   /* Skip next three lines what are comments. */
      for(m = 0; m < 3; m++) while(fgetc(rptr) != '\n');

   /* Load data at current time-step. */
      for(i = 0; i < g->nr[isec]+2*NrG; i++) for(j = 0; j < g->nth+2*NthG; j++) for(k = 0; k < g->nz[isec]+2*NzG; k++)
          {
           n = 0; while((str[n] = fgetc(rptr)) != '\n') n++; str[n] = '\0';

           m = 4;
#          ifdef SOLVE_ENERGY
           m = 5;
#          endif

           s = parse_str(str, '\t', 1, m, &n);

           if(n != m)
            {
             fprintf(g->lptr, "ERROR:: In 'init_cond()':  Reading data at current time, t = %E sec: the recovery file '%s' is not in proper format.\n", g->t, recovf);
             exit(-1);
            }

           g->u[isec][i][j][k] = str2num(s[0], 3); free(s[0]);
           g->v[isec][i][j][k] = str2num(s[1], 3); free(s[1]);
           g->w[isec][i][j][k] = str2num(s[2], 3); free(s[2]);
           g->p[isec][i][j][k] = str2num(s[3], 3); free(s[3]);

#          ifdef SOLVE_ENERGY
           g->T[isec][i][j][k] = str2num(s[4], 3); free(s[4]);
#          endif

           free(s);
          }
     }


/* Load data at outlet-boudary of the SECTION-2 for PREVIOUS time-steps. These are used
   for the applciation of Orlansky BC to determine phase velocities at outlet. */
   for(tn = 1; tn < tn_BACKUP; tn++)
     {
     /* tn = 0: t = t; It is filled later in "orlansky_boundary()" using current time-step data.
        tn = 1: t = t - dt;
        tn = 2: t = t - 2*dt; etc. */

   // Skip further three lines what are comments.
      for(m = 0; m < 3; m++) while(fgetc(rptr) != '\n');

      for(i = 0; i < g->nr[2]; i++) for(j = 0; j < g->nth; j++) for(k = 0; k < zCells_BACKUP; k++)
       {
        n = 0; while((str[n] = fgetc(rptr)) != '\n') n++; str[n] = '\0';

        m = 3;
#       ifdef SOLVE_ENERGY
        m = 4;
#       endif

        s = parse_str(str, '\t', 1, m, &n);

        if(n != m)
            {
             fprintf(g->lptr, "ERROR:: In 'init_cond()':  Reading data at past time-instants at outlet: the recovery file '%s' is not in proper format.\n", recovf);
             exit(-1);
            }

        g->outBC.u[tn][i][j][k] = str2num(s[0], 3); free(s[0]);
        g->outBC.v[tn][i][j][k] = str2num(s[1], 3); free(s[1]);
        g->outBC.w[tn][i][j][k] = str2num(s[2], 3); free(s[2]);
#       ifdef SOLVE_ENERGY
        g->outBC.T[tn][i][j][k] = str2num(s[3], 3); free(s[3]);
#       endif

        free(s);
       }
     }

   free(str);

   if(fgetc(rptr) != EOF)
       {
        fprintf(g->lptr, "ERROR:: In 'init_cond()': The recovery file, %s, doesn't seem to be in correct format.\n", recovf);
        exit(-1);
       }

   fclose(rptr);

#endif

   return;
  }
